function [fbasedk] = flowbased3dk(k,dx,dy,dz,transmi)

% permute(k,[2 3 1]);
rot90_3D(k,1,1);

%Number of cells in each direction 
[nc_j,nc_i,nc_k] = size(k);

%Test if need upscaling 
if numel(k) == 1
    fbasedk = k;
    return
end

if transmi == 1 % for transmissibility upscaling
    dx = 1;
    dy = 1;
    dz = 1;
end

%Contour i = 1 correction
if nc_i == 1
    k = cat(2,k,k);
    [nc_j,nc_i,nc_k] = size(k);
    dx = dx/2;
end

%Number of points in each direction
np_i = nc_i-1;
np_j = nc_j+1;
np_k = nc_k+1;

%Total points
np = np_i * np_j * np_k;

p1 = 1; %  %Boundary condition;
% pn = 0; %  %Boundary condition; %no use

%finite difference matrix
%prelocating nodal equations vectors
top = zeros(np,1);

back = zeros(np,1);
left = zeros(np,1);
center = zeros(np,1);
right = zeros(np,1);
front = zeros(np,1);

buttom = zeros(np,1);

%prelocating left term vector
b = zeros(np,1);

%%
for x = 1:np
    
    [ci,cj,ck] = ponto(x,np_i,np_j);
    
    %     %%
    if ck == 1
        
        %%
        if cj == 1
            
            ktop = 0;
            kback = k(cj,ci,ck);
            kleft = 0;
            
            kright = mean([k(cj,ci,ck) k(cj,ci+1,ck)]);
            kfront = k(cj,ci+1,ck);
            kbuttom = mean([k(cj,ci,ck) k(cj,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            %%
        elseif cj == np_j
            
            ktop = 0;
            kback = k(cj-1,ci,ck);
            kleft = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck)]);
            
            kright = 0;
            kfront = k(cj-1,ci+1,ck);
            kbuttom = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            
            %%
        else
            %%
            
            ktop = 0;
            kback = mean([k(cj-1,ci,ck) k(cj,ci,ck)]);
            kleft = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck)]);
            
            kright = mean([k(cj,ci,ck) k(cj,ci+1,ck)]);
            kfront = mean([k(cj-1,ci+1,ck) k(cj,ci+1,ck)]);
            kbuttom = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck) ...
                k(cj,ci,ck) k(cj,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            
        end
        
        %%
    elseif ck == np_k
        
        %%
        if cj == 1
            
            ktop = mean([k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kback = k(cj,ci,ck-1);
            kleft = 0;
            
            kright = mean([k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kfront = k(cj,ci+1,ck-1);
            kbuttom = 0;
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            %%
        elseif cj == np_j
            
            ktop = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            kback = k(cj-1,ci,ck-1);
            kleft = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            
            kright = 0;
            kfront = k(cj-1,ci+1,ck-1);
            kbuttom = 0;
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            
            %%    %         elseif (ci == 1 || ci == np_i) && (cj ~= 1 || cj ~= np_j)
        else
            %%
            ktop = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1) ...
                k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kback = mean([k(cj-1,ci,ck-1) k(cj,ci,ck-1)]);
            kleft = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            
            kright = mean([k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kfront = mean([k(cj-1,ci+1,ck-1) k(cj,ci+1,ck-1)]);
            kbuttom = 0;
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            
        end
        %%
    else
        
        %%        
        if cj == 1;
            
            ktop = mean([k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kback = mean([k(cj,ci,ck) k(cj,ci,ck-1)]);
            kleft = 0;
            
            kright = mean([k(cj,ci,ck) k(cj,ci+1,ck) ...
                k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kfront = mean([k(cj,ci+1,ck) k(cj,ci+1,ck-1)]);
            kbuttom = mean([k(cj,ci,ck) k(cj,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
         %%   
        elseif cj == np_j
            %%
            ktop = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            kback = mean([k(cj-1,ci,ck) k(cj-1,ci,ck-1)]);
            kleft = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck) ...
                k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            
            kright = 0;
            kfront = mean([k(cj-1,ci+1,ck) k(cj-1,ci+1,ck-1)]);
            kbuttom = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
           
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end
            
        else
            %%
            ktop = mean([k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1) ...
                k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kback = mean([k(cj-1,ci,ck) k(cj,ci,ck) ...
                k(cj-1,ci,ck-1) k(cj,ci,ck-1)]);
            kleft = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck) ...
                k(cj-1,ci,ck-1) k(cj-1,ci+1,ck-1)]);
            
            kright = mean([k(cj,ci,ck) k(cj,ci+1,ck) ...
                k(cj,ci,ck-1) k(cj,ci+1,ck-1)]);
            kfront = mean([k(cj-1,ci+1,ck) k(cj,ci+1,ck) ...
                k(cj-1,ci+1,ck-1) k(cj,ci+1,ck-1)]);
            kbuttom = mean([k(cj-1,ci,ck) k(cj-1,ci+1,ck) ...
                k(cj,ci,ck) k(cj,ci+1,ck)]);
            
            [top(x),back(x),left(x),center(x),right(x),front(x),buttom(x)] = ...
                vetor(dx,dy,dz,ktop,kback,kleft,kright,kfront,kbuttom);
            
            if ci == 1
                
                b(x) = -back(x);
                back(x) = 0;
                
            end
            if ci == np_i

                front(x) = 0;
                
            end          
            
        end
        %%
        %%
    end
end
%%
% Sparse Matrix Formed From Diagonals
M = spdiags([buttom front right center left back top],...
    [-np_i*np_j -np_j -1 0 1 np_j np_i*np_j],np,np);

% spy(M)

%%   
%Solving system by gauss-elimination 
p = M\b;

%Reshaping to 3D matrix
% p = reshape(p,[np_j np_i np_k]);  %there is no need

%Add - Boundry Condition
% p = [p1*ones(np_j,1,np_k) p pn*ones(np_j,1,np_k)];  %there is no need  
%%
%Flux calculate - we need to get the flux passing by to determine fbasedk
flux = zeros(np,1);

for z = 0:np_k-1
    
   zcont = z*(np_i*np_j)+1 : z*(np_i*np_j)+np_j;
   
    flux(zcont) = b(zcont).*(p(zcont)-p1);
           
end
%%
totalflux = sum(flux);

if transmi == 0
    %permeability    
    fbasedk = totalflux*(nc_i*dx)/(nc_j*dy*nc_k*dz);
else 
    %transmissiblity ===> T = q :P
    fbasedk = totalflux;
end

%%
end